home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / tutorial < prev    next >
Text File  |  1994-05-19  |  8KB  |  321 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26. /* tutorial.c 10/12/1993 */
  27.  
  28. /* routines from Chapter 1 of Meschach */
  29.  
  30. static char rcsid[] = "$Id: tutorial.c,v 1.3 1994/01/16 22:53:09 des Exp $";
  31.  
  32. #include <math.h>
  33. #include "matrix.h"
  34.  
  35. /* rk4 -- 4th order Runge--Kutta method */
  36. double rk4(f,t,x,h)
  37. double t, h;
  38. VEC    *(*f)(), *x;
  39. {
  40.    static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL;
  41.    static VEC *temp=VNULL;
  42.    
  43.    /* do not work with NULL initial vector */
  44.    if ( x == VNULL )
  45.      error(E_NULL,"rk4");
  46.  
  47.    /* ensure that v1, ..., v4, temp are of the correct size */
  48.    v1   = v_resize(v1,x->dim);
  49.    v2   = v_resize(v2,x->dim);
  50.    v3   = v_resize(v3,x->dim);
  51.    v4   = v_resize(v4,x->dim);
  52.    temp = v_resize(temp,x->dim);
  53.  
  54.    /* register workspace variables */
  55.    MEM_STAT_REG(v1,TYPE_VEC);
  56.    MEM_STAT_REG(v2,TYPE_VEC);
  57.    MEM_STAT_REG(v3,TYPE_VEC);
  58.    MEM_STAT_REG(v4,TYPE_VEC);
  59.    MEM_STAT_REG(temp,TYPE_VEC);
  60.    /* end of memory allocation */
  61.  
  62.    (*f)(t,x,v1); /* most compilers allow: "f(t,x,v1);" */
  63.    v_mltadd(x,v1,0.5*h,temp);    /* temp = x+.5*h*v1 */
  64.    (*f)(t+0.5*h,temp,v2);
  65.    v_mltadd(x,v2,0.5*h,temp);    /* temp = x+.5*h*v2 */
  66.    (*f)(t+0.5*h,temp,v3);
  67.    v_mltadd(x,v3,h,temp);        /* temp = x+h*v3 */
  68.    (*f)(t+h,temp,v4);
  69.    
  70.    /* now add: v1+2*v2+2*v3+v4 */
  71.    v_copy(v1,temp);              /* temp = v1 */
  72.    v_mltadd(temp,v2,2.0,temp);   /* temp = v1+2*v2 */
  73.    v_mltadd(temp,v3,2.0,temp);   /* temp = v1+2*v2+2*v3 */
  74.    v_add(temp,v4,temp);          /* temp = v1+2*v2+2*v3+v4 */
  75.    
  76.    /* adjust x */
  77.    v_mltadd(x,temp,h/6.0,x);     /* x = x+(h/6)*temp */
  78.    
  79.    return t+h;                   /* return the new time */
  80. }
  81.  
  82.  
  83.  
  84. /* rk4 -- 4th order Runge-Kutta method */
  85. /* another variant */
  86. double rk4_var(f,t,x,h)
  87. double t, h;
  88. VEC    *(*f)(), *x;
  89. {
  90.    static VEC *v1, *v2, *v3, *v4, *temp;
  91.    
  92.    /* do not work with NULL initial vector */
  93.    if ( x == VNULL )        error(E_NULL,"rk4");
  94.    
  95.    /* ensure that v1, ..., v4, temp are of the correct size */
  96.    v_resize_vars(x->dim, &v1, &v2, &v3, &v4, &temp, NULL);
  97.  
  98.    /* register workspace variables */
  99.    mem_stat_reg_vars(0, TYPE_VEC, &v1, &v2, &v3, &v4, &temp, NULL);
  100.    /* end of memory allocation */
  101.  
  102.    (*f)(t,x,v1);             v_mltadd(x,v1,0.5*h,temp);
  103.    (*f)(t+0.5*h,temp,v2);    v_mltadd(x,v2,0.5*h,temp);
  104.    (*f)(t+0.5*h,temp,v3);    v_mltadd(x,v3,h,temp);
  105.    (*f)(t+h,temp,v4);
  106.    
  107.    /* now add: temp = v1+2*v2+2*v3+v4 */
  108.    v_linlist(temp, v1, 1.0, v2, 2.0, v3, 2.0, v4, 1.0, VNULL);
  109.    /* adjust x */
  110.    v_mltadd(x,temp,h/6.0,x);     /* x = x+(h/6)*temp */
  111.    
  112.    return t+h;                   /* return the new time */
  113. }
  114.  
  115.  
  116. /* f -- right-hand side of ODE solver */
  117. VEC    *f(t,x,out)
  118. VEC    *x, *out;
  119. double    t;
  120. {
  121.    if ( x == VNULL || out == VNULL )
  122.      error(E_NULL,"f");
  123.    if ( x->dim != 2 || out->dim != 2 )
  124.      error(E_SIZES,"f");
  125.    
  126.    out->ve[0] = x->ve[1];
  127.    out->ve[1] = - x->ve[0];
  128.    
  129.    return out;
  130. }
  131.  
  132.  
  133. void tutor_rk4()
  134. {
  135.    VEC        *x;
  136.    VEC        *f();
  137.    double     h, t, t_fin;
  138.    double     rk4();
  139.    
  140.    input("Input initial time: ","%lf",&t);
  141.    input("Input final time: ",  "%lf",&t_fin);
  142.    x = v_get(2);    /* this is the size needed by f() */
  143.    prompter("Input initial state:\n");    x = v_input(VNULL);
  144.    input("Input step size: ",   "%lf",&h);
  145.    
  146.    printf("# At time %g, the state is\n",t);
  147.    v_output(x);
  148.    while (t < t_fin)
  149.    {
  150.       /* you can use t = rk4_var(f,t,x,min(h,t_fin-t)); */
  151.       t = rk4(f,t,x,min(h,t_fin-t));   /* new t is returned */
  152.       printf("# At time %g, the state is\n",t);
  153.       v_output(x);
  154.    }
  155. }
  156.  
  157.  
  158.  
  159.  
  160. #include "matrix2.h"
  161.  
  162. void tutor_ls()
  163. {
  164.    MAT *A, *QR;
  165.    VEC *b, *x, *diag;
  166.    
  167.    /* read in A matrix */
  168.    printf("Input A matrix:\n");
  169.    
  170.    A = m_input(MNULL);     /* A has whatever size is input */
  171.    
  172.    if ( A->m < A->n )
  173.    {
  174.       printf("Need m >= n to obtain least squares fit\n");
  175.       exit(0);
  176.    }
  177.    printf("# A =\n");       m_output(A);
  178.    diag = v_get(A->m);
  179.    /* QR is to be the QR factorisation of A */
  180.    QR = m_copy(A,MNULL);
  181.    QRfactor(QR,diag);   
  182.    /* read in b vector */
  183.    printf("Input b vector:\n");
  184.    b = v_get(A->m);
  185.    b = v_input(b);
  186.    printf("# b =\n");       v_output(b);
  187.    
  188.    /* solve for x */
  189.    x = QRsolve(QR,diag,b,VNULL);
  190.    printf("Vector of best fit parameters is\n");
  191.    v_output(x);
  192.    /* ... and work out norm of errors... */
  193.    printf("||A*x-b|| = %g\n",
  194.       v_norm2(v_sub(mv_mlt(A,x,VNULL),b,VNULL)));
  195. }
  196.  
  197.  
  198. #include "iter.h"
  199.  
  200.  
  201. #define N 50
  202. #define VEC2MAT(v,m)  vm_move((v),0,(m),0,0,N,N);
  203.  
  204. #define PI 3.141592653589793116
  205. #define index(i,j) (N*((i)-1)+(j)-1)
  206.  
  207. /* right hand side function (for generating b) */
  208. double f1(x,y)
  209. double x,y;
  210. {
  211.   /* return 2.0*PI*PI*sin(PI*x)*sin(PI*y); */
  212.    return exp(x*y);
  213. }
  214.  
  215. /* discrete laplacian */
  216. SPMAT *laplacian(A)
  217. SPMAT *A;
  218. {
  219.    Real h;
  220.    int i,j;
  221.    
  222.    if (!A)
  223.      A = sp_get(N*N,N*N,5);
  224.  
  225.    for ( i = 1; i <= N; i++ )
  226.      for ( j = 1; j <= N; j++ )
  227.      {
  228.         if ( i < N )
  229.       sp_set_val(A,index(i,j),index(i+1,j),-1.0);
  230.         if ( i > 1 )
  231.       sp_set_val(A,index(i,j),index(i-1,j),-1.0);
  232.         if ( j < N )
  233.       sp_set_val(A,index(i,j),index(i,j+1),-1.0);
  234.         if ( j > 1 )
  235.       sp_set_val(A,index(i,j),index(i,j-1),-1.0);
  236.         sp_set_val(A,index(i,j),index(i,j),4.0);
  237.      }
  238.    return A;
  239. }
  240.  
  241. /* generating right hand side */
  242. VEC *rhs_lap(b)
  243. VEC *b;
  244. {
  245.    Real h,h2,x,y;
  246.    int i,j;
  247.    
  248.    if (!b)
  249.      b = v_get(N*N);
  250.  
  251.    h = 1.0/(N+1);      /* for a unit square */
  252.    h2 = h*h;
  253.    x = 0.0;
  254.    for ( i = 1; i <= N; i++ ) {
  255.       x += h;
  256.       y = 0.0;
  257.      for ( j = 1; j <= N; j++ ) {
  258.     y += h;
  259.     b->ve[index(i,j)] = h2*f1(x,y);
  260.      }
  261.    }
  262.    return b;
  263. }
  264.    
  265. void tut_lap()
  266. {
  267.    SPMAT *A, *LLT;
  268.    VEC *b, *out, *x;
  269.    MAT *B;
  270.    int num_steps;
  271.    FILE *fp;
  272.  
  273.    A = sp_get(N*N,N*N,5);
  274.    b = v_get(N*N);
  275.  
  276.    laplacian(A);
  277.    LLT = sp_copy(A);
  278.    spICHfactor(LLT);
  279.  
  280.    out = v_get(A->m);
  281.    x = v_get(A->m);
  282.  
  283.    rhs_lap(b);   /* new rhs */
  284.    iter_spcg(A,LLT,b,1e-6,out,1000,&num_steps);
  285.    printf("Number of iterations = %d\n",num_steps);
  286.  
  287.    /* save b as a MATLAB matrix */
  288.  
  289.    fp = fopen("laplace.mat","w");  /* b will be saved in laplace.mat */
  290.    if (fp == NULL) {
  291.       printf("Cannot open %s\n","laplace.mat");
  292.       exit(1);
  293.    }
  294.    
  295.    /* b must be transformed to a matrix */
  296.    
  297.    B = m_get(N,N);
  298.    VEC2MAT(out,B);
  299.    m_save(fp,B,"sol");  /* sol is an internal name in MATLAB */
  300.  
  301. }
  302.  
  303.  
  304. void main()
  305. {
  306.    int i;
  307.  
  308.    input("Choose the problem (1=Runge-Kutta, 2=least squares,3=laplace): ",
  309.      "%d",&i);
  310.    switch (i) {
  311.     case 1: tutor_rk4(); break;
  312.     case 2: tutor_ls(); break;
  313.     case 3: tut_lap(); break;
  314.     default: 
  315.       printf(" Wrong value of i (only 1, 2 or 3)\n\n");
  316.       break;
  317.    }
  318.  
  319. }
  320.  
  321.